1. Field of the Invention
The invention relates to the field of control engineering. It relates to a method and device for estimating a state of a system.
2. Discussion of the Background
Tracking the temporal performance of a technical system, for example a production plant, a power plant or an electromechanical device, requires continuous determination of system variables, for example temperatures, pressures, mass flows, outputs. With the aid of such variables and with the aid of values of control inputs, the system can be controlled and the future performance of the system can be predicted, optimized and held within safe limits.
The measurement of all the required system variables is, however, neither possible nor expedient in technical or economic terms. Consequently, specific system variables have to be estimated with the aid of measurements of other variables. All estimates are affected by random errors because of measuring noise, variations in the sensors and non-observed disturbances. Tracking and predicting system performance, and regulating and controlling the system are based on these faulty estimates.
System variables are generally estimated using a Kalman filter. If the system has linear dynamics, the filter determines from the measurements a point in the space of the system variables. However, since the physical laws acting in the system, for example the laws of conservation of mass and energy, are not all linear, they are not obeyed by the estimate of the filter. Consequently, the estimates are not of adequate quality for accurate monitoring.
For this reason, methods exist or coordinating measured data, also termed data reconciliation methods, which ensure that at least the measured data obey prescribed physical laws or satisfy boundary conditions.
The term xe2x80x9cstatexe2x80x9d is not used below in the narrower sense of dynamic system theory, that is to say not only for the minimum information for describing the temporal system performance, but in the wider sense for all system variables which are of interest.
Thus, the aim is to determine a value x of a d-dimensional state vector X of a system. The information contained in this state is usually redundant because of physical relationships in the system. This redundancy corresponds to a set of boundary conditions corresponding to the c-dimensional equation
g(x)=0.xe2x80x83xe2x80x83(1)
Let a measurement Y of the state X be disturbed by an instance of independent additive noise W. The latter represents sensor noise, for example. These relationships are reproduced by FIG. 1. Thus, it holds that
Y=X+W.xe2x80x83xe2x80x83(2)
Here and below, upper-case letters W, X, Y denote random vectors, that is to say vectors of random variables, while lower-case letters w, x, y denote values which these random vectors assume.
Because of the boundary conditions (1), the prior probability density of X satisfies
px(x)xe2x88x9dxcex4(g(x)),
where xcex4(xc2x7) is the Dirac distribution. Furthermore, because X and W have been assumed to be independent, the likelihood function is
py|x(y|x)=pw(yxe2x88x92x),xe2x80x83xe2x80x83(3)
where Pw(xc2x7) is the probability density of the measuring noise. Given a measurement, in the a posteriori probability density of X in accordance with Bayes theorem in which the random variable Y assumes a value y, that is to say that Y=y, is given by                                           P                          X              |              Y                                ⁡                      (                          x              |              y                        )                          =                                                                              p                  x                                ⁡                                  (                  x                  )                                            ⁢                                                p                  w                                ⁡                                  (                                      y                    -                    x                                    )                                                                                    ∫                                  -                  ∞                                ∞                            ⁢                                                                    p                    x                                    ⁡                                      (                    x                    )                                                  ⁢                                                      p                    w                                    ⁡                                      (                                          y                      -                      x                                        )                                                  ⁢                                  xe2x80x83                                ⁢                                  ⅆ                  x                                                              .                                    (        4        )            
Given that Y=y, the maximum a posteriori (MAP) estimate, that is to say the most probable value of X is therefore                                                                                                               x                    ^                                    MAP                                ⁡                                  (                  y                  )                                            ⁢                              =                ^                            ⁢                            ⁢                                                                    arg                    ⁢                                          xe2x80x83                                        ⁢                    max                                                        {                                                                  x                        |                                                  g                          ⁡                                                      (                            x                            )                                                                                              =                      0                                        }                                                  ⁢                                                      p                    x                                    ⁡                                      (                    x                    )                                                  ⁢                                                      p                    w                                    ⁡                                      (                                          y                      -                      x                                        )                                                                                                                          =                            ⁢                                                                    arg                    ⁢                                          xe2x80x83                                        ⁢                    min                                                        {                                                                  x                        |                                                  g                          ⁡                                                      (                            x                            )                                                                                              =                      0                                        }                                                  -                                  log                  ⁢                                      xe2x80x83                                    ⁢                                                            p                      x                                        ⁡                                          (                      x                      )                                                                      -                                  log                  ⁢                                      xe2x80x83                                    ⁢                                                            p                      w                                        ⁡                                          (                                              y                        -                        x                                            )                                                                                                                              (        5        )            
For example, let the measurement be modeled by additive Gaussian noise with a mean value of zero and non-correlated components, that is to say
Wxcx9cN(0;xcexa3w),xe2x80x83xe2x80x83(6)
The covariance matrix being equal to
xcexa3w=diag("sgr"12, . . . , "sgr"m2)xe2x80x83xe2x80x83(7)
Assuming a constant prior density (uninformative prior) px(xc2x7) for X, (5) yields the maximum likelihood (ML) estimate, {circumflex over (x)}ML(y) as solution of                                           min            x                    ⁢                                                    (                                  y                  -                  x                                )                            T                        ⁢                                          ∑                W                                  -                  1                                            ⁢                              xe2x80x83                            ⁢                              (                                  y                  -                  x                                )                                                    =                              min            x                    ⁢                                    ∑                              i                =                1                            m                        ⁢                          xe2x80x83                        ⁢                                                            (                                                            y                      i                                        -                                          x                      i                                                        )                                2                            /                              σ                i                2                                                                        (        8        )            
under the boundary condition that
g(x)=0.xe2x80x83xe2x80x83(9)
This corresponds to the classical data reconciliation method such as is known, for example, from the article xe2x80x9cData Reconciliation, Progress and Challengexe2x80x9d; Cameron M. Crowe; J. Proc. Cont. Vol. 6, No. 2/3, pp. 89-98, 1996; Elsevier Science Ltd. The quality of the method is inadequate for certain applications.
It is therefore the object of the invention to create a method and a device of the type mentioned at the beginning which remove the above-mentioned disadvantages.
In the method according to the invention, after a measurement a state of a technical system is estimated with the aid of new measured data from the system, from a measure of the accuracy of the measured data and from physical boundary conditions. In this case, a preferably immediately preceding state, a measure of its accuracy and a dynamic model of the system are also taken into account, and a measure of the accuracy of the newly estimated state is determined.
The method according to the invention renders it possible to incorporate not only the static physical boundary conditions, but also their dynamic development. Consequently, the accuracy of the estimated state and the quality of regulation or control based thereon are enhanced. For example, the robustness of an estimate with respect to noise-induced and other outliers is substantially enhanced. Such an outlier would cause a false estimated value {circumflex over (x)}ML(y) in the case of conventional methods.
A system is to be understood as a technical process of arbitrary type, for example a plant used in production engineering, an electromechanical device or a power plant, in particular a thermal power plant for generating electric energy.
Here and below the term xe2x80x9cstatexe2x80x9d is not used in the narrower sense of dynamic system theory, that is to say it is used not only for the minimum information for describing the temporal system performance, but in the wider sense for all system variables which are of interest. For the sake of simplicity, the term xe2x80x9caccuracyxe2x80x9d is used at some junctures instead of the term xe2x80x9cmeasure of accuracyxe2x80x9d.
In a preferred variant of the method according to the invention, the accuracy of the system variables is used to determine their uncertainty, for example in the form of probability intervals, and account is taken of this in regulating the system. Consequently, expected operating costs can be determined and optimized continuously and comprehensively. In particular, operating risks can be reduced, for example by taking account of the uncertainty of a critical system variable when this variable must be kept in a specific range.
A device according to the invention has a means for estimating state variables under boundary conditions, a means for predicting state variables, a means for estimating accuracies and a means for estimating an accuracy of predictions.
A preferred embodiment of the device according to the invention further has means for determining uncertainties and means for optimizing and regulating the performance of the system.